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ABSTRACT 


This report is an improvement to ISAP-1, "SRB Vorticity- 
Acoustic Coupled Instability Analysis", September, 1986. 

Included in this report are the automatic generation of all input 
data for grid configuration, boundary conditions for coupled 
acoustic and vortical field calculations, transformation of all 
dimensions to a parametric form, resulting in flexibility for the 
user to define the size of the problem (geometric configurations) 
with reduction in storage (15-65%) and computer time (50-75%) . 
Additional research is required for the following areas: 

(1) effects of turbulence, (2) nonlinear wave oscillations, and 
(3) chemistry upon combustion instability. 
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I . INTRODUCTION 


This report represents an improved version of "SRB 
Vorticity-Acoustic Coupled Instability Analysis - ISAP-l, 
September, 1986. The many basic changes to the original code 
include the automatic generation of all the input data for grid 
structure, boundary conditions, and coupling between the flow 
field and the acoustic/ vortical field. In addition, all the 
dimensions in the program were transformed into a parametric 
form. These new parameters will enable the user to control the 
computer memory storage and the program execution time by 
specifying different sets of parameters and different geometric 
configurations. Also presented is the comparison between the 
results of the original program and those of the new one, along 
with recommendations to the user and additional research 
requirements . 

As noted in the earlier report (ISAP-l) , unstable waves may 
occur as a result of acoustic and/or vortical (hydrodynamic) 
oscillations. If these two different types of waves are coupled 
together, their physical interactions lead to extremely 
complicated phenomena. Theoretically, there exists an infinite 
number of frequencies for both acoustic and vortical 
oscillations. Realistically, however, only a limited number of 
combined frequencies are excited. Our objective is to determine 
the combined nature of acoustic and vortical frequencies at which 
instabilities may arise. This subject is important in rocket 
motor combustion chambers when the vortical field is coupled with 
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acoustic pressure oscillations. In the past, the acoustic 
combustion instability was studied independently of the vortical 
instability induced by vortex motions. This report is intended 
to combine the two different sources of energy everywhere within 
the spatial domain and to determine the effect of one upon the 
other. This can be achieved by calculating the mean flow 
velocities and vorticities and their fluctuating parts of 
velocities and vortices, as well as the fluctuating pressure. 

To elucidate this coupling mechanism, the acoustic wave 
equation and the perturbed vortical transport equation are 
solved, being combined with the results of mean flow calculations 
from the Navier-Stokes system by means of finite elements. With 
these data, growth constants are calculated and stability 
boundaries determined. Contributions to stability and/or 
instability from various sources such as combustion, convection 
(flow turning) , and viscous damping on propellant surfaces and 
energy convection, momentum convection, momentum viscous damping, 
and dissipative energy from the interior domain are separately 
identified. It is also found that stability boundaries for 
coupled acoustic and vortical oscillations are somewhat similar 
to the classical hydrodynamic stability boundaries, but they 
occur in the form of multiple islands. 
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II. MESH GENERATION 

In the original version of the program, the input data was 
read externally, and this procedure required extensive 
preparation for any changes in the input constants and/or the 
configuration of the field. This process is complex and time 
consuming. Consequently, a mesh generation routine is added to 
the program. 

The program has two major parts: (1) the flow field 
calculations and (2) the acoustic and vortical calculations from 
which the stability integrals and growth constants are derived. 
Therefore, a complete set of input data is required for each 
part. 


A. Flow Field 

The flow field calculations include velocities and pressure 
with grid configurations coarser than those in ISAP-1. Figure la 
shows the original grid in ISAP-1. It is apparent that the 
element sizes are smaller near boundaries (4), (5), and (6) in 

Fig. 2, but these reductions were set arbitrarily. The mesh 
generation routine has a much better approach. It reduces the 
grid size logarithmically near the same boundaries. This will 
enhance the flexibility of the program. The generated arrays are 
the following: 

(a) NENN - the element connectivity matrix; it sets the global 
nodal values of the nodes of each element. 

(b) XX, YY - the coordinates of each global node. 
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(c) NU, NV, and ND - the global nodes for the Dirichlet boundary 
condition adjustment for the U, V velocities and pressure, 
respectively. 

(d) UB, VB, and PB - the boundary values associated with NU, NV, 
and ND, respectively. 

(e) UU, W, PP - the velocities in the x-direction, y-direction, 
and the pressure, respectively. 

B. Acoustic £Qd Vortic al Fi e ld 

This is the smaller grid shown in Fig. 2b, d. The two grids 
are interconnected through the ICON matrix. Again this 
connectivity was previously set arbitrarily. While in the mesh 
generation routine, the new grid is set through ICON, such that 
the distances between the nodes at each boundary are as equal as 
possible. The arrays for this field are: 

(a) NENL - the element connectivity matrix for the smaller grid. 

(b) ICON - the interconnectivity matrix between the flow field 
and the acoustic/vortical field. 

(c) NODE - the adjustment matrix for the boundary conditions in 
the computation of vortical modes. 

(d) FFNX, FFNY - the direction cosines at the boundaries in the 
x and y directions, respectively. 

(e) NBQ - the number of elements at each boundary. 

(f) ABN, ANN - the admittance at the burning surface and nozzle, 
respectively . 

(g) NC - the connectivity matrix of the boundary elements. 
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C. Initial flaw Constants 

For each problem , there are some basic constants that are 
needed to define it. They deal with the geometry, flow 
properties, error, time-step, and convergence parameters. Many 
of these constants, such as REN, GAMMA, DT, ERROR, ITMAX, and 
NPT, have already been defined in the original report. In order 
to enhance the flexibility of the program, a few additional 
constants dealing with the geometry of the problem were utilized 
in the mesh generation routine. They are the following: 

(a) XMIN, XMAX, YMIN, YMAX - the limiting values of the domain. 

(b) XO, YO - the coordinates of the corner node (see Fig. 1). 
These geometric values must be redefined by the user every 

time the configuration of the problem changes. 
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III. PARAMETRIC DIMENSIONS 


* The original program was bounded by a set of fixed 
parameters which controls the memory storage requirements. For 
the flow field, the program used 832 grid points and 767 
elements; for the acoustic/vortical field, it used 40 grid points 
and 27 elements. 

In order to make the programs more flexible, parametric 
dimensions were added. This addition will allow the user to 
choose the memory storage needed for different executions. This 
procedure requires the definition of the parameters and the 
respective arrays. 

A. Parameters 

A total of 10 parameters is needed to maximize the 
flexibility of the program. All but two of these parameters are 
geometric; they will generate both grids. Of the remaining two 
parameters, one is fixed and the other is floating (not fixed) . 
(a) Geometric Parameters : 

Four parameters are needed for each grid. 

1. Flow field : 

The four parameters, shown in Fig. 2, are defined as 
follows: 

NELF1 - the number of boundary elements above the corner node 
in the y-direction. 

• the number of boundary elements below the corner node 
in the y-direction. 


NELF2 
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NETP1 - the number of boundary elements left of the corner node 
In the x-direction. 

NETP2 - the number of boundary elements right of the corner 
node in the x-direction. 

2. Acoustic/vortical field: 

The four parameters associated with this field are similar 
to those of the flow field with their respective positions. 

NALF1, NALF2, NAT PI , and NATP2 are the number of elements on 
boundaries (5), (1), (6), and (4), respectively (see Fig. 2). 

(b) Floating Par amete r: 

The parameter MI is affected by other parameters and 
constants and is a direct result of the calculations. In the 
original code, the value of MI was set at 65. However, this 
value can be much lower (see Table 3). 

From the program, the needed value of MI is augmented for 
each St (Strouhal number) > .01. Those values of St < .01 will 
produce oscillations which are negligible relative to the 
acoustic oscillations and, therefore, are eliminated to avoid 
unnecessary calculations. As the Reynolds number increases, the 
vortical oscillations expand further in the field, resulting in 
higher values of St at additional nodes; therefore, the needed 
value of MI will increase. 

However, in order to limit the memory storage requirements, 
MI must be minimized. This procedure will be achieved by 
changing MI relative to the changes in REN and the geometric 
parameters (see Section V for the recommended values of MI) . 
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(c) Fixed Ear a nte ter : 

The parameter NBP, which Is always equal to 6, represents 
the* 6 boundaries of the domain . It must be noted that the 
remaining parameters are derived directly from the 10 preset 
parameters already defined. 

B. Arrays 

The size of all the different arrays can be set either by a 
DIMENSION statement or by a COMMON block. When the parameters 
were added to the program, two other steps had to be taken: 

(1) All the COMMON blocks were eliminated from the code. This 
step required a change in the way the subroutines are 
called, by including all the variable arrays that were 
otherwise included in the COMMON statements. 

(2) All the arrays must be referenced in the MAIN program. This 
step was required because the dimensions of all the arrays 
are generated in the PARAMETER statements in the MAIN 


program. 
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IV. COMPARISON OF RESULTS 

* The computer memory storage and the execution time are the 
important outcome of this study. Therefore, to compare these 
effects, the results will be displayed in two ways: (1) by 
comparing the original code and the revised code and (2) by 
comparing the original grids' results with those of the reduced 
ones. 

First of all, the basic assumptions and flow constants are 
presented: 

1. The Reynolds number, REN - 10 3 , 10 4 , 10 s for different runs. 

2. The specific heat ratio, GAMMA - 1.2. 

3. The time step for the calculation of the mean flow field, 

DT - 1. 

4. The convergence error for mean flow, ERROR - .001. 

5. The maximum number of iterations for flow field 
calculations, ITMAX - 30. 

6. The domain of calculations: XMIN ■ 0, XMAX » 10, YMIN - 0, 
YMAX - 1.5. 

7. The corner node coordinates: X0 ■ 2, YO - 1. 

8. Boundary and initial conditions: With reference to Fig. 1, 

burning takes place only on boundary (6) ; therefore, an 
instantaneous flux normal to this boundary appears at time 
T * 0. Then, the velocity in the negative y-direction is 
set at a dimensionless value of v - -.01, only at the nodes 
of boundary (6), and excluding the end points of this 
boundary. All other flow variables are set equal to zero 



everywhere in the domain. (The value of -.01 was taken from 
the original code and is a function of the burning rate at 
the surface of the solid fuel) . 

9. The normal vectors at the boundaries: At boundaries (1) - 
(6), the values of the normal vectors are FFNX = 1., 0., l. f 

0. , 1., 0.; FFNY =0., 1., 0., 1., 0., 1., respectively. 

10. Admittance: Only at the burning surface (boundary 6), ABN = 

1. ; otherwise, it is equal to zero. Only at the nozzle 
(boundary 3), AAN =1.; otherwise, it is equal to zero. 

These assumptions and flow constants were kept fixed (except 

for REN) throughout the calculations. Now, we can proceed to the 
comparison of the results. 

A. Original Grids 

It is clear that with the original grids (ISAP-1) , the 
storage required to generate the data in the revised code will 
exceed the storage needed to read the data from an outside 
source. Therefore, the total memory requirements for the revised 
code exceed that of the original code by about 2%. 

Concerning the run time requirements, on the other hand, the 
revised code used about 2% less time than that of the original 
code, even though it went through 5 additional values of Strouhal 
numbers, i.e., vortical calculations (see Table 1). 

NOTE: The value of MI has been defined earlier, but this value 

is also dependent on the systems on which these programs are 
running. The original code gave a value of MI = 15 on the IBM 
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3084 and a value of MI - 16 on the UNIVAC. This phenomenon is a 
result of the eigenvalue solver routine. This routine is very 
sensitive to any small variations in its input data. To 
eliminate the sensitivity problem, these codes should be 
transformed to operate in double-precision, while reducing the 
ERROR value and increasing ITMAX. These new additions will 
increase the storage requirements and run time by over 50%. 

The flow field results, as shown in Fig. 3, prove that the 
revised code presents the same exact flow field as that of the 
original code. In the figure, the original code results are 
displayed in the upper half of each section for (a) the pressure, 
(b) the velocity in the x-direction, and (c) the velocity in the 
y-direction. These results were expected because of the 
linearity of the solution. However, the difference between the 
values of MI was not expected to be as large. 

Even though the revised version of the program made it more 
flexible and easy to operate, it still has not made any major 
contribution to reduce the memory requirements nor the necessary 
run time. Therefore, we will present the other aspects of 
possible improvements which reduce the memory storage and run 
time. But, as we will see, they will slightly affect the 
computational results. 

B. Reduced Grids 

As a part of the flexibility of the revised code, the user 
has the option of changing the number of nodes and elements used 
in the program. Figure 1 shows the different grids used for 
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comparison purposes. The number of nodes and elements are 
tabulated in Table 2. As seen in Table 2, if the number of 
elements at the boundaries is reduced by 1/4, the total number of 
nodes and elements is reduced by about 1/2. As a result of the 
reduced grids, the required memory storage and run time will 
definitely be reduced. Then, there is a need to find what 
percentage reduction is achieved and how it will affect the 
results . 

Table 3 shows the different reductions achieved for required 
memory and run time at different Reynolds numbers. The storage 
requirement is reduced by 14% and 66% with coarser grids for thm 
flow field and the flow and A/V fields, respectively, while the 
execution time is reduced between 56% and 76%. From these 
values, we can conclude that the flow field calculations affect 
mostly the run time required, while the acoustic/vortical (A/V) 
field affects the storage size of the program. On the other 
hand, the computational results of such reductions are slightly 
affected. Figures 4-6 show the flow field variables for the 
reduced grid (the lower part of each section) , as comapred to the 
original grid, for Re - 10 3 , 10 4 , and 10 5 , respectively. The 
difference between the results becomes more obvious as Re gets 
larger. Consequently, at high values of the Reynolds number 
(Re), the coarse mesh does not produce accurate results. 
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V. RECOMMENDATIONS 

* Two typos of recommendations are necessary to conclude this 
report: (1) those directed towards the users of the program and 
(2) those needed for a more realistic modeling of the problem at 
hand. 

A. To the user: 

(1) Grids: This program does not contain a turbulent flow model 

to accommodate for the higher values of Re. Also, there is 
a limit, in reducing the grids, at which the results would 
become insignificant. Therefore, a coarse mesh, such as (6, 
9, 13, 20), may be used at Re < 10 3 ? but, at higher values 
of Re, a more refined mesh is needed. Recommended values 
for the parameters: NELF1 > 6, NELF2 > 9, NETP1 > 13, 

NETP2 > 20, 2 < NALF1, NALF2, NATP1 , NATP2 <4. It is 
imperative that the last 4 parameters be > 2; otherwise, the 
calculations will become senseless. 

(2) The value of MI: As stated earlier, the needed value of MI 

is calculated within the program. The maximum allowable 
value of MI is equal to the number of A/V modes: MI m#x » NT. 
In order to control the memory storage requirements, MI can 
be specified smaller than NT. As shown in Table 3, the 
needed values of MI increase as the Reynolds number 
increases . 

Therefore, it is recommended that the user adjust the value 
of MI, as needed, to correspond with the values of the Reynolds 
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number and the appropriate geometric parameters. The recommended 
values of HI are: MI (ne , d#d) < MI < NT, where MI (n ,, ded) are the 

values listed in Table 3. 

B. Additional Recommendations : 

There are many areas where this program can be improved, 
some of which can be done in the near future, others will be the 
topics of further study and research. 

(1) Double precision, smaller DT, smaller ERROR, and larger 
ITMAX: these extra steps are needed for accuracy and to 
avoid the sensitivity in the eigenvalue solver. 

(2) Turbulent model: this model will be used for the relatively 
high values of Re. 

(3) Higher order approximations of the pressure and the 
vorticity. 

(4) The study of the effect of compressibility. 

(5) The addition of chemistry to the model by adding the species 
and energy equations. 
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Fig 2 . Parameters, boundaries, and domain of calculations. 




(C) 


Fif 3. Comparison between the flow field results of the 
original code (upper part) and those of the revised 
code (lower oart) , for (a) the pressure, (b) the 
U-velocity, and (c) the V-velocity. (Re=1000.) 



(C) 


Fig 4. Comparison of the flow field variables between the 

original mesh (upper part) and the reduced mesh (lower 
part) : (a) the pressure, (b) the U-velocity, (c) the 
V-velocity; Re= 1000. 
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(c) 


Fig 5. Comparison of the flow field variables between the 

original mesh (upper part) and the reduced mesh (lower 
part) : (a) the pressure, (b) the U- velocity, (c) the 
V-velocity; Re= 10000. 



(c) 


Fig 6. Comparison of the flow field variables between the 

original mesh (upper part) and the reduced mesh (lower 
part) : (a) the pressure, (b) the U-velocity, (c) the 
V-velocity; Re* 100000. 


A/V Field Flow Field 
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Table 1. Comparison between the storage and runtime of the 
original code and the revised one, at Re= 10^ 
(original grids) 



Memory 

Requirement 

Bytes 

Fraction 
of the 
original 

Run-time 
(UN I VAC) 
min 

Fraction 
of the 
original 

MI 

(needed) 

ORIGINAL 

3023004 

1 

226 

1 

16 

REVISED 

3109452 

1.02 

223 

.98 

21 


Table 2. The reduced grids. 


Parameters 


No. of 
Nodes 


No. of 
Elements 


Original 


(8,13,17,26) 


832 


767 


Reduced 


(6,9,13,20) 


466 


417 


Original 


( 2 , 3 , 4 , 3 ) 


40 


27 


Reduced 


( 2 , 2 , 2 , 2 ) 


21 
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Flow § A/V Flow field Original 


Table 3. Storage and run-tiim 
grids and the reduc 
an initial assumpti 


Re 

Memory 
Requi remen 

Bytes 

10 3 

3109452 

10 4 

3109452 

10 5 

3109452 

10 3 

2682808 

10 4 

2682808 

10 s 

2682808 

10 3 

1077816 

4 

10 4 

1077816 

10 5 

1077816 


Fractio 


comparison between the original 
ones. (MI was set at 65 as 
). 



65 


.270 


14 
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Listing of the Revised Code. 
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